Graphite is composed of many graphene monolayers, where carbon atoms in the layer are arranged into a hexagonal lattice connected by strong covalent bonds and only weak Van der Waals (vdW) force exists between layers. The study of phonon cross-plane (CP) thermal conduction in graphite is helpful in revealing the phonon transport mechanism in vdW layered material in the vertical direction^[27] or across the vdW interface between individual nanostructures or nanostructure and substrate/host material.^[2]

It is generally accepted that graphite has a very low CP thermal conductivity. The measurements^[3,4] performed decades ago already showed that the graphite CP thermal conductivity is only about 7 W/mK. However, there are still many debates about the mechanism of phonon CP transfer. For example, the average values of mean free path (MFP) in CP direction are significantly different in different models. In the classic kinetic theory based on an isotropic assumption, is less than one nanometer in Ref. [5], and about 6.7 nm in Ref. [6]. As a comparison, Sadeghi et al.^[7] performed a calculation by using the full phonon dispersion of graphite based on the density functional perturbation theory, and estimated at 20 nm under 300 K. In the molecular dynamics simulation, is on the order of several hundred nanometers.^[8] Harb et al.^[9] measured a CP thermal conductivity of a graphite thin film to be 35 nm by time resolved x-ray diffraction, and reported a value of about 0.7 W/mK much lower than the value of bulk graphite, which suggests that MFP in CP direction can be longer than 35 nm. Fu et al.^[10] measured the CP thermal conductivity of thin graphite flakes with different thickness values by using a differential three-omega method, and then estimated at around 200 nm at room temperature. In a similar experiment carried out by Zhang et al.,^[11] through using the time-domain thermoreflectance technique, is found to be a few hundred nanometers at room temperature and a much narrower distribution than in isotropic crystals. It should be noted that MFP cannot be measured directly but it can be estimated by the MFP reconstruction method to link the measured thermal conductivities to the MFP spectrum.^[12,13] However, this is a challenging problem because different MFP spectra may lead to the same thermal conductivity.

According to the experimental measurements, Yang et al.^[14] suggested that the high frequency phonons can propagate only in the basal plane but not in the CP direction because of the weak vdW interaction between basal planes. Therefore, these high frequency phonons should not be included in the calculation, and, as a result, could be over 100 nm at room temperature. Recently, a similar idea is found in the work of Han et al.^[15] on the assumption that high frequency phonons are confined in the graphite layer, and have to exchange energy with the low-frequency phonons under 4 THz which are the only phonons allowed to transfer in the CP direction. Furthermore, taking into account only a small part of volumetric heat capacity contributed from these low frequency phonons, they predict a c-axis MFP of 165 nm at 300 K based on a modified one-dimensional kinetic theory.

There are still many controversies and debates about the contributions from different phonon modes. In this paper, we will give a wave vector and polarization/branch resolved analysis of all the phonon modes in Brillouin zone (BZ) directly, with the help of the almaBTE first-principles phonon solver.^[16] The rest of this paper is organized as follows. In Section 2, the method used in the calculation is briefly given. In Section 3, we calculate the phonon dispersion, thermal conductivity and specific heat at constant volume, including a comparison with measurements to make sure that we have a good starting point. The CP thermal conductivity, specific heat and MFP of different branches in different part of BZ are then analyzed. Finally, some conclusions are drawn from the present study.

Under the relaxation time approximation (RTA) and in the framework of Boltzmann transport equation, the CP thermal conductivity can be expressed as^[17]

The integrations in Eqs. (1), (2), and (4) run over the entire BZ of bulk graphite crystal, which is hexagonal, as shown in Fig. 1. In the almaBTE package,^[16] the BZ is discretized into a regular wave vector grids, denoted by green dots in Fig. 1. There are grids with denoting the in-plane (IP) direction and the CP direction. Each grid contains a small volume fraction of the BZ, and has a wave vector , frequency for the phonon polarization p according to the phonon dispersion relation. With the almaBTE package, we can furthermore show the distributions of many physical quantities, such as relaxation time , specific heat , thermal conductivity, etc. in the wave vector grids. Hence the integration over BZ is transformed into summations over all phonon polarizations and wave vector grids. The wave vector grid of 36×36×12 as illustrated in Fig. 1 is used to fulfill the convergence requirements of and κ.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Schematic illustration of BZ of graphite with wave vector grids denoted by green dots. The black dots refer to the irreducible wave vector grids. Ga and Gc denote length of reciprocal-lattice vector along the plane and CP direction, respectively. R is radius of a cylinder with axis coinciding with Γ–A. The cylinder is used to choose phonon modes in the center of BZ.

The calculated graphite phonon dispersion is denoted by dots in Fig. 2. The triangles, circles and squares refer to the experiment results.^[19–21] These phonons are in the high-symmetry directions Γ K, Γ M, and Γ A in BZ as shown in Fig. 1. The calculated results accord well with the measurements, which presents a solid foundation for further study.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Phonon dispersion in graphite along high-symmetry direction. Blue dots are for calculations. Triangles and circles are for inelastic x-ray scattering measurements from Refs. [19] and [20] respectively. Squares are for neutron scattering data from Ref. [21]. The solid lines are used to show the two branches ZA and , which are discussed in detail in this paper.

The building block of graphite is the graphene monolayer, which has two atoms per unit cell. These layers are bound together by a weak vdW force. There are six phonon branches for graphene. Each of these six branches splits into two branches for graphite which has four atoms per unit cell.^[22] Therefore, graphite has twelve branches, i.e., three acoustic (A) branches (LA, TA, ZA) and nine optical (O) branches. Except for ZA and , the remaining branches are double degenerate almost in the entire BZ. Herein L, T, and Z represent longitudinal polarization, in-plane and out-of-plane transversal polarization, respectively. Strictly speaking, these labels hold true only along the high-symmetry directions in BZ, such as Γ K and Γ M, and are misleading for arbitrary wave vectors. However, these labels are useful in classifying the branches. For the ZA banch and branch denoted by lines in Fig. 2, the atoms in the graphite layers vibrate along the out-of-plane direction. These two branches are split from graphene ZA (flexural) branch, which also has out-of-plane vibrations and exhibits a quadratic dispersion near Γ^[22] However, as a three-dimensional material, graphite has wave vectors which are not confined in the in-plane directions. Therefore, ZA and are, in fact, respectively longitudinal acoustic and optical branch for phonons with wave vectors along Γ A in BZ^[23], as shown in Fig. 2(b).

Based on the phonon dispersion over the entire BZ grids in Fig. 1, the calculated CP thermal conductivity κ_total represented by the line with square symbols slightly overestimates the measurement κ_experiment^[24] denoted by triangles in Fig. 3(a) in a temperature range of 100 K–1000 K. This result is reasonable because the only scattering taken into account is the three-phonon process.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (a) Comparison between CP thermal conductivity measurements (triangles)[24] and the calculated results of κtotal (line with squares), including its ZA and branch contributions. (b) Temperature-dependent different phonon percentage contributions.

Although there are 12 branches, only two branches, i.e., ZA and represented by solid lines in Fig. 2, have significant contributions to CP thermal conduction, as reported by Paulatto et al.^[22] The temperature dependence of (line with dots) and (line with diamonds) are illustrated in Fig. 3(a), and their percentage contributions are shown in Fig. 3(b). The sum of the partial contributions due to ZA and phonons denoted by the line with circle is almost 100% of total thermal conductivity, and slightly decreases as temperature rises because high temperature can excite more high-frequency phonons belonging to other branches to participate in thermal conduction. From 200 K to 1000 K, about 80% and 20% of contributions come from ZA and phonons, respectively.

Based on the phonon dispersion over the entire BZ, the specific heat at constant volume of graphite is calculated as illustrated in Fig. 4. The C_total (line with dots) takes into account all the branches and accords well with the experiment values (line with triangles).^[25] The ZA and branch, denote respectively by the line with x symbols and the line with circles almost have the same specific heat in a temperature range of 100 K–1000 K. The dotted lines identify the high temperature limits for different phonon branches in which each phonon mode simply contributions to the specific heat the same constant value as required by the Dulong and Petit law.^[26] As illustrated in Fig. 2, the ZA and branch each have a low maximum phonon frequency. Therefore, as temperature rises, it is easy to satisfy which will lead to the specific heat per mode in Eq. (2) . At 100 K, the ZA branch and branch have the dominant contributions to the specific heat. A further rise in the temperature will excite more high phonons belonging to other branches, which leads to a dramatic increase in the total specific heat C_total, while C_ZA and begin to approach to the high temperature limits and cannot increase significantly. This situation means that the sum of C_ZA and represented by the dashed line in Fig. 4 can only contribute a small part of the total specific heat as temperature rises.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Plots of specific heat at constant volume of graphite versus temperature. The line with dots takes into account all the branches. ZA branch and branch and their sum are represented by line with x symbols, line with circles, and dashed line, respectively. The line with triangles is for experiment results from Ref. [25]. Dotted lines represent high-temperature limits.

The results shown in Figs. 3 and 4 reveal that in the view of thermal conduction in the CP direction, except for ZA and , the remaining phonon branches occupy most of thermal energy but are trapped in graphite layers above 300 K. Polarized in the CP direction, ZA and phonons act as the dominant interlayer thermal energy carriers. Graphite, as a vdW layered material, has a weak interaction between layers. The vdW force is sensitive to the distance between atoms. Atoms in the neighbor layers have the smallest distance with relatively strong vdW interaction, and vibrate in the CP direction to change the interatomic distance dramatically in the ZA mode and mode. This kind of vibration can benefit the interlayer transfer of thermal energy.

The wave vector resolved plots of CP thermal conductivity of ZA branch and branch are illustrated in Fig. 5. The results are for the irreducible wave vector grids highlighted as black dots in Fig. 1. These grids represent the discrete irreducible BZ, the smallest wedge reduced from BZ by all of the symmetries. Each grid occupies a small volume of BZ, in which phonons have contributions to CP thermal conductivity, represented by the color of dots in Fig. 5. The exact value of the contribution can change with the density of grids, or the grid volume. However, Figure 5 reveals the distribution of contributions from phonons in different parts of BZ.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Distributions of CP thermal conductivity contribution from (a) branch ZA and (b) branch over the irreducible wave vector grids at 300 K.

For ZA branch, Figure 5(a) reveals that the phonons in the grids along the Γ–A line have major contributions, and significantly larger than in the other parts of BZ. The CP thermal conductivity of phonons has a significantly smaller value than that of ZA phonons, and is relatively evenly distributed in the identical grids as shown in Fig. 5(b). Moreover, the phonons in the vicinity of Γ–A line and near to point A are more likely to participate in the CP thermal conduction. It should be noted that because only a few grids have very large contributions for ZA phonons, is only about four times the value of at 300 K as shown in Fig. 3(b).

To find which parts of phonons in BZ are critical for CP thermal conduction, a cylinder shape with its axis coinciding with Γ–A is defined in Fig. 1. The relation between the normalized cumulative CP thermal conductivity at 300 K and the cylinder radius R is illustrated in Fig. 6. The lines with circles, triangles and squares represent the case in which κ takes into account the contribution from all branch phonons, ZA phonons, and phonons in the cylinder, respectively. When , phonons in the cylinder transfer about 80% thermal energy in the CP direction. A similar conclusion was obtained by Paulatto et al.^[22] In addition, ZA phonons have the dominant contribution, especially when R is small. It appears that ZA phonons are the only thermal energy carriers along the Γ–A line, and contributes about 30% of the CP transport. The phonons have significantly less contribution than ZA phonons, and the contribution is approximately proportional to the cylinder radius until it reaches its maximum value 20% at . Phonons near the edge of BZ ( ) have almost no contribution.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Plots of normalized cumulative CP thermal conductivity at 300 K versus radius R of cylinder in BZ with axis coinciding with Γ–A.

Graphite is famous for its anisotropic properties. Therefore, it is reasonable to explore the direction dependence of phonon transport. Herein, we divide the entire BZ into six sectors according to the wave vector angle (or the group velocity angle with respect to the graphite layer normal, and assess the contribution due to phonons in these sectors. Figure 7 shows the slices of the iso-energy surfaces in the A–Γ–M plane which is divided uniformly into six sectors by the dash-dotted lines with respect to . The vertex angle of each sector at the point Γ is 15°. The arrows indicate the directions of the group velocities which are normal to the iso-energy surfaces. The length of arrows represents the velocity magnitude. There are six velocity arrow groups identified by different colors according to the angle specified in the legend below the Figure.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Slices of the iso-energy surfaces of (a) ZA and (b) phonon in A–Γ–M plane. The plane is divided into 6 sections by dash dotted lines. according to angle between wave vector and Γ–A direction. Group velocities denoted by arrows are divided into 6 groups identified by different colors with legend shown below the Figure according to angle between velocity direction and Γ–A direction.

As an acoustic phonon, the ZA mode has a zero frequency at the Γ point, which is the center of an ellipsoid iso-energy surface in the low frequency range as shown in Fig. 7(a). Furthermore, the ellipsoid surface grows in the graphite flat hexagonal prism shape BZ as frequency rises, and transforms into a cylinder surface coaxial with Γ–A. As shown in Fig. 7(b), the branch has similar iso-energy surfaces, but its lowest frequency of optical phonons is not zero and is at the A point.

The anisotropic ellipsoid shape and the cylinder shape of iso-energy surface are attributed to the strong bonding between atoms in the graphite layer accompanied with the weak vdW interaction between layers. The group velocity arrows normal to the iso-energy surfaces in Fig. 7 indicate the direction of energy flow. The anisotropic shape of iso-energy surface results in the existence of regions that each have a small curvature, where phonons tend to transfer energy in nearly the same direction, instead of along the wave vector direction. This phenomenon, especially when phonon ballistic transport plays a remarkable role at low temperature, is called phonon focusing, which means that the uniform angular distribution of wave vectors gives rise to a nonuniform distribution of group velocity vectors in an anisotropic medium and can be experimentally observed.^[27]

Because of the anisotropic shapes of iso-energy surfaces, the sectors are not the same as the sectors in Fig. 7. As a result, Figure 8 shows that there is significant difference in the contribution of CP thermal conductivity between the phonons in sectors (dashed line with solid symbols) and those in sectors (solid line with hollow symbols) at 300 K. The research about the wave vector angle dependence shows that the relatively large contributions come from two sectors: [0°,15°) and [60°,75°), and the contribution of phonons increases with increasing until a maximum value only about 8% in sector [60°,75°) is reached.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. CP thermal conductivity contributions at 300 K in and sectors defined in Fig. 7.

Figure 8 reveals that phonons with group velocity close to CP direction and IP direction play an important role in CP thermal conduction. When the velocity vector angle , the red arrows are distributed in a large part of BZ in the vicinity of Γ–A line, in which the largest wave vector angle is about 45° for the ZA branch in Fig. 7(a), and about 60° for the branch in Fig. 7(b). These phonons denoted by red arrows are apparently focused in the CP direction. When all branches are taken into account, the contribution (solid line with hollow dots in Fig. 8) due to phonons with has almost 44%, including 40% from the ZA phonons. If phonons with are also included, the total contribution will rise to 50%. By comparison, phonons with wave vector only transfer about 33% thermal energy in the CP direction, indicated by dashed line with solid dots in Fig. 8.

Phonons with blue arrows in Fig. 7 have group velocity vector angle . A rough estimation tends to underestimate the contributions from these phonons. However, these blue arrows do have a component in the CP direction and the more important fact is that these phonons occupy at least 3/4 of BZ where phonon modes are distributed uniformly. When all the branches are considered, the contribution in is about 23%, including 14% from ZA phonons and 8% from phonons.

Based on alamaBTE,^[16] results can be given for each grid which represents a small part of BZ shown in Fig. 1. Therefore, the MFP projected in the CP direction for each mode can be calculated. The largest value of of ZA phonon and phonon are found to be about and 52 nm at 300 K. It is apparent that the transport properties of phonons can change dramatically and the MFP can span several orders of magnitude.^[8] To assess MFP in a big picture view, the average MFP defined in Eq. (4) is calculated, and the value is about 15 nm at 300 K, which is comparable to 20 nm estimated by Sadeghi et al. by using a theoretical method based on the full phonon dispersion.^[7] The values of are about 23 nm and 8 nm respectively for branch ZA and branch . It is apparent that the average MFP is far lower than the maximum mode MFP . The reason is that as illustrated in Figs. 3 and 4, the ZA branch and branch have the major contributions to CP thermal conduction, but have a small part of specific heat. Even for one branch, each different phonon mode has very significantly different . As a result, the average value calculated based on Eq. (4) is small.

Figure 9 shows the plots of angle-dependent at 300 K averaged in different sectors defined in Fig. 7. As discussed previously, phonons have small contributions to CP thermal conduction. One reason revealed by Fig. 9 is that phonons have less than 20 nm in all the sector (dashed line with solid squares) or sectors (solid line with hollow squares). The dominant thermal energy carriers, ZA phonons that have long are gathered around the Γ–A line. Therefore, the value of all branches or ZA branch alone approximately tends to decrease as or increases. When the wave vector angle , of ZA branch denoted by dashed line with solid triangles has a large value of 737 nm. If all the branches are considered, in sector decreases to 375 nm (dashed line with solid dots). When the group velocity angle , the of all branches (solid line with hollow circles) and ZA branch alone (solid line with hollow triangles) are 76 nm and 112 nm, respectively. As the longitudinal acoustic phonons, these ZA phonons with long focus in the CP direction, their contribution to thermal conductivity is about 40% in as discussed previously in Fig. 8.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Plots of average projected MFP at 300 K in CP direction in different sections defined according to angle ( or range.

In summary, we have investigated the phonon transport along the CP direction of graphite from first-principles calculations and in the single-mode relaxation time approximation. The calculated phonon dispersion, CP thermal conductivity and specific heat for graphite are in good agreement with measurements. In a temperature range from 200 K to 1000 K, about 80% and 20% of contributions to CP heat transfer come respectively from ZA and phonon branch, which are polarized in the CP direction, and can be regarded as longitudinal acoustic and optical branch for phonons with wave vectors along the Γ–A line in BZ, respectively. Except for ZA and ZO’, the remaining phonon branches occupy most of thermal energy, but almost cannot transfer between graphite layers above 300 K. The phonons in the central 16% of BZ contribute 80% of CP transport. Phonons with group velocity angle with respect to Γ–A less than 30° have 50% of CP transport at 300 K. The largest MFP in CP direction at 300 K is found to be about and 52 nm for ZA and phonon, respectively. The phonons with long MFP are gathered around Γ–A. When the group velocity angle , the average MFP in CP direction at 300 K for ZA branch is about 112 nm.

The author thanks Li Deyu for helpful discussion.